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The bispectrum of single-held inhat ionary trajectories in which the speed of sound of the in¬ 
hationary trajectories c s is constant but not equal to the speed of light c = 1 is explored. The 
trajectories are generated as random realisations of the Hubble Slow-Roll (HSR) hierarchy and the 
bispectra are calculated using numerical techniques that extend the work of Jl]. This method allows 
for out-of-slow-roll models with non-trivial time dependence and arbitrarily low c s . The ensembles 
obtained using this method yield distributions for the shape and scale-dependence of the bispec¬ 
trum and their relations with the standard inhationary parameters such as scalar spectral tilt n s 
and tensor-to-scalar ratio r. The distributions demonstrate the squeezed-limit consistency relations 
for arbitrary single-held inhat ionary models. 


I. INTRODUCTION 

Current observations of the universe suggest that its 
density perturbations, to a good approximation, can be 
considered as a realisation of a correlated Gaussian statis¬ 
tic and are very close to but not exactly scale independent 
HE]. This scale dependence is characterised by the mea¬ 
surement of the scalar spectral index n s = 0.968 db 0.006 
[4] which agrees well with the framework of the early 
universe undergoing a phase of quasi-de Sitter expansion 
that resulted in correlated, super-horizon scaled curva¬ 
ture perturbations to the background metric. The stan¬ 
dard, and the most commonly accepted, explanation for 
both the origin of the perturbations and the reason for 
the quasi-de Sitter expansion is the presence of a scalar 
held known as the inhaton whose potential energy dom¬ 
inates the Hubble equation and whose spatial fluctua¬ 
tions seed the curvature perturbations that later drive 
all structure formation EHm. 

One of the main issues facing efforts aimed at under¬ 
standing the nature and origin of the inhaton is that 
many classes of different inhationary models predict ob¬ 
servables such as n s and r that are in broad agreement 
with observations (see for example [T8H23] ). With the 
hnal analysis of Planck data imminent and the combined 
Planck-BICEPII/Keck analysis [24] conhrming that r 
was in fact not detected in the BICEPII data [25] this 
situation may become the status quo for the foreseeable 
future. This will be the case unless tensor modes, in the 
form of r 0, are detected by the next generation of sub¬ 
orbital Cosmic Microwave Background (CMB) experi¬ 
ments, or, non-Gaussianity is measured. In the former 
case, discernment between different inhationary models 
may also require the measurement of the spectral tilt 
of tensor modes n t which is challenging due to the cos¬ 
mic variance effect on the largest scales where the tensor 
mode signal is clearest. 

A detection of non-Gaussianity, in the form of a non¬ 
zero bispectrum [23 [27] or un-connected contributions 
to higher order moments, may then provide the key to 
uncovering the origin of the inhaton. Non-Gaussianity is 
necessarily present in the universe since general relativity 
is a non-linear theory and even if the inflation were driven 


by a single, free, scalar held it would still interact with 
gravity giving rise to a non-zero bispectrum. In general, 
the non-Gaussianity of less standard models of inflation, 
particularly ones that predict low tensor contributions 
with r —>> 0, tends to be large and potentially measurable 
in the near future. 

The bispectrum is the third-order moment of the cur¬ 
vature perturbation in Fourier space and is expected to 
be the easiest non-Gaussian signal to measure as it is 
both the lowest order component in the perturbation and 
has no Gaussian counterpart. Observational bounds are 
often quoted in terms of the scale-free amplitude /nl 
[28] , a dimensionless quantity which is typically of order 
the Slow-Roll (SR) parameter e ^ 10 -2 for simple inha¬ 
tionary models [29., 30]. For more complicated models, 
it is possible to generate a larger /nl while maintain¬ 
ing n s « 1 and much effort has been spent constructing 
such models in the hope that a large non-Gaussianity is 
detected (see [2TH211 E6l EU (32] for some examples). 

Within the context of single held models, there are 
a couple of possibilities. One is to break the slow- 
roll approximation temporarily by introducing a feature 
[33[[34], such as a bump, in the inhaton potential E (</>). 
A second is to use a non-canonical kinetic term for the 
scalar held [22J |3T[ [32j [35] . This involves adding extra 
derivatives d^cj) as interactions for the held. One phys¬ 
ical consequence of this is that the scalar perturbations 
typically propagate at a new sound speed c s < 1 and it 
is these models that will be considering in this work. 

In this work, for simplicity, we restrict ourselves to 
the case of a constant c s =/= 1, reserving arbitrary time- 
dependent sound speeds for future work. We calculate 
the bispectrum of these models numerically, allowing for 
high values of c~ 2 — 1 and combine this with a Monte 
Carlo approach for sampling inhationary models. We 
analyse in detail the exact scale and shape dependence of 
such models, verifying our results by demonstrating the 
squeezed-limit consistency relation for very small sound 
speeds and large SR parameters. 

This paper is organised as follows; In Section [II] we 
summarise the framework and parameters required for 
the calculation of the bispectrum and briefly discuss the 
Monte Carlo generation of inhationary trajectories using 
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the Hamilton-Ja cobi fo rmalism discussed in more detail 
in pQ. In Section III A we give an overview of the numer¬ 
ical calculation of the power spectrum before proceed¬ 
ing to the calculation of the bispectrum in Section |IIIB| 
We summarise our results and consistency checks in Sec¬ 
tion IV before finally concluding in Section [V| 


II. MONTE-CARLO APPROACH TO 
SAMPLING TRAJECTORIES 


This Hamilton-Jacobi (HJ) formalism j36ll39] . and its 
role in numerical inflation was discussed at length in | 1 ] 
and we refer the reader to that work for an extended 
discussion. Here we summarise the method. In the HJ 
formalism the dynamics of an inflating cosmology can be 
captured entirely by considering the Hubble parameter, 
iJ( 0 ) as a function of the inflaton field value </> and by 
considering a hierarchy of Hubble Slow-Roll (HSR) pa¬ 
rameters defining the hierarchy if derivatives of H with 
respect to <fi. 

We extend this formalism by introducing an arbitrary, 
but constant sound speed c s ^ 1. Following mMM 
we consider actions of the form 


J d 4 Xy/=gC, 

(i) 

M vl 


~X~R + P(x, 4>) i 

(2) 


(3) 


where M p i is the Planck mass, R is the Ricci scalar, and 
is the inverse space-time metric. The Lagrangian 
density C in the action above describes a perfect fluid 
with pressure P(X, </>) and energy density p = 2 XPx — P 
where Px = dP/dX. The speed of sound, c s , is defined 
as 


2 _ Px_ _ Px 
Px Px + 2 XP xx 


(4) 


For constant c s this can be treated as a differential equa¬ 
tion for P(X, (j)). Using the initial condition P(X, 0) = 
X — V (</>) when c s = 1 one obtains 


P(X^) = 


2 c 2 


1 


W 2 


- V(ct >). 


(5) 


The equation of motion for <p differs from the canonical 
case so the original definitions of the HSR parameters in 
the HJ formalism should be altered accordingly. How¬ 
ever, one can still define e —foldings A, the Hubble rate 
H(t) and its time derivatives independently of the dy¬ 
namics of the inflation. That is 


a(N) = 

e\ 

( 6 ) 

H(N) = 

d d N 

(7) 

a d t 

e(N) = 

din H 

( 8 ) 

d N ’ 


where a is the scale factor and overdots denote differen¬ 
tiation with respect to cosmic time t. The HSR param¬ 
eters can now be defined so that they correspond to the 
HJ formalism HSR parameters in the limit where c s = 1 

aI\ 

— = [le + (l -l)n\ l \- i+ 1 A, (9) 

where 1 A = 77, 2 A = £. 

The values of ^A at the end of inflation at N = N tot 
can be drawn randomly to sample the distribution of con¬ 
sistent inflationary trajectories as described in [1]. The 
sound speed will not affect the time dependence of these 
parameters so it will not play an explicit role in the sam¬ 
pling of trajectories . In practice the random sampling is 
achieved by drawing the following set of parameters with 
uniform distributions (flat prior) in the intervals 

l X=[-l,l]xe~ sl (10) 

N tot = [60,80]+In A, (11) 


where / > 0. In addition since we draw samples at the end 
of inflation we fix the value of the l = 0 HSR parameter 
°A = e(N tot ) = 1 . 

In (10), x and s are parameters that specify the scaling 
of the uniform prior range with l and can be used to inves¬ 
tigate the dependence of our final results on the assumed 
priors. The random sampling of N tot represents the un¬ 
certainty in the total duration of the post-inflationary 
reheating phase and the constant A is related to the nor¬ 
malisation of H which will be discussed shortly. For¬ 
mally one would need to evolve an infinite number of 
l X parameters to sample the space of all possible H(N) 
functions. In practice this is not possible and one must 
truncate the series at some finite order L max . We define 
Lmax such that P max HSR parameters includes e(A) e.g. 
Umax = 2 corresponds to e(A) and 77 (A) with all other 
^A = 0 identically. Once random values of l X have been 
drawn the entire inflationary trajectory can be obtained 
by integrating the background equations of motion suffi¬ 
ciently far back in the past to cover the required number 
of e- foldings given by N tot . 


III. COMPUTATIONAL METHOD 

The calculation of the bispectrum relies on the same 
basic building blocks as the calculation of the primor¬ 
dial power spectrum. In addition the bispectrum is often 
compared to the spectral tilt of the power spectrum and 
the squeezed limit consistency condition is a valuable tool 
for checking the numerical method. We therefore give a 
brief review the calculation of the power spectrum as the 
first step in the numerical calculation of the bispectrum. 


A. Computation of the power spectrum 

We choose a gauge where the inflaton perturbation 
J 0 (t, x) = 0 and the spatial metric is given by gij = 
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FIG. 1: Dependence of /nl on the damping factor S when n = 1 for squeezed (left) and folded (right) configurations. 
For this trajectory c~ 2 = 3. The red-solid and blue-dashed lines show = 10 -5 (Mpc) -1 and 
kmax = 10 -2 (Mpc) -1 respectively. For large 5 the damping factor is too large affecting the horizon crossing 
behaviour and the oscillations provide no contribution, producing a smooth curve. For small S the oscillations are 
not sufficiently suppressed producing noise. Ideally /nl should converge with decreasing S in some sense to its true 
value before the noise begins to dominate. There is an indication of this in the right panel at S ~ 0.1. Unfortunately 
for the squeezed limit, the amplitude of /nl is too small relative to the noise to extract any reasonable result. To 
make matters worse, depending on the shape the optimum S changes by an order of magnitude. Also note noise 

begins at larger S for larger k. 


a 2 (t)e 2( ^ t,x ) Sij . This defines the comoving curvature per¬ 
turbation £(£, x). The primordial power spectrum of the 
curvature perturbation is then 

(C kl (* k2 ) = (27r) 3 <5( 3 )(k 1 +k 2 )P c (fci), (12) 

where k is the wavevector of the Fourier mode and k = 
|k|. These modes satisfy the Mukhanov-Sasaki equation 
hzieb which, with our choice of variables becomes 


d 2 4 

dTV 2 


+ (3 +1 


0 , 

lr,) AN + 


is' 


(13) 


To obtain the power spectrum we simply require the 
freeze-out value of (k when the mode crosses the sound- 
horizon, i.e. 


m o = ici 2 


c s k<^aH 


(14) 


wavenumber k is set by requiring that c s k = AaH(N *.) 
where A 1. In practice this means that the integra¬ 
tion is started at successively later times as k increases. 
This avoids unnecessary computational steps at smaller 
scales. 


Assuming H(N ) varies slowly enough, each mode will 
evolve for roughly In A e-folds before they cross the 
sound-horizon and freeze out. The earliest mode of in¬ 
terest to freeze out will be k m i n so we choose A^ min = 0, 
i.e. TV = 0 is defined such that c s k m i n = AaH(N = 0) 
and we then apply (15) to this mode. This means the 
k m in mode will cross the sound-horizon at N c « In A and 
we can then use the standard analytical result relating 
H to the amplitude of the power spectrum to normalise 
H . In practice, during the backwards integration of the 
HSR parameters, we apply a normalisation condition on 
H such that 


notice that for theories where the speed of sound and light 
are not equivalent the horizon set by the speed of sound 
is the relevant scale beyond which freeze-out occurs. 

We apply the usual Bunch-Davies initial conditions 
[42] when the mode is deep inside the sound-horizon 


(k 


2 M p ia V ke 


—ic s kr 


(15) 


where r is conformal time defined through dN/dr = aH. 

We impose initial conditions © at different e-folds for 
each mode k. This ensures all modes are sufficiently deep 
inside the sound-horizon at the start of the forward inte¬ 


gration of (13). The starting e- folds, 7V&, for mode with 


H{N C ) = 2ir^2^(NcjA s M pl , (16) 


where A s is conventional the normalisation of the dimen¬ 
sionless primordial curvature power spectrum. In the 
usual power law convention for the form of the power 
spectrum A s is employed as 


k 3 P c (k) = a 2 s 


n s — 1 


(17) 


A similar procedure can be carried out for the calcu¬ 
lation of the gravitational wave spectrum which is un¬ 
affected by c s . The analogues of (13) and (15) are are 
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FIG. 2: Dependence of /nl on the damping factor S when n = 3 for squeezed (left) and folded (right) configurations. 
For this trajectory c~ 2 = 3. The red-solid and blue-dashed lines show = 10 -5 (Mpc) -1 and 
kmax — 10 -2 (Mpc) -1 respectively. By choosing n > 1 the suppression is weighted more towards the early time 
oscillations and less on the horizon crossing time. Practically this pushes the noise back to very small values of S 
allowing /nl to converge to its true value. The acceptable range of S is also much wider solving the shape 
dependence problem. One could choose n 1 but most of the time c s k/aH > 1. Therefore to prevent damping at 
horizon crossing 8 must be reduced to compensate and this method can only be pushed so far. In practice we found 
n = 3 to be sufficient. Note for large k the noise still arises at larger S. 


identical to the standard case with c s 


d 2 h k dh k f k 

dtp +{3 - e hN + {^s 

hk 


hk = 0 , 

1 e~ ikr 
M p i aV2k 


(18) 

(19) 


A complication that arises due to the sound and light 
horizon not being the same is that scalar and tensor 
modes freeze out at different times so one must be sure 
that the Bunch-Davies conditions are applied when both 
modes are sufficiently deep inside their respective hori¬ 
zons. In principle the power spectrum must converge in 
the limit A oo therefore the answer should not depend 
on whether the Bunch-Davies conditions are applied ear¬ 
lier to one mode with respect to another as long as both 
modes are sufficiently deep inside their respective hori¬ 
zons. In practice this means nothing needs to be changed. 
If c s k = AaH then we know k > AaH as c s < 1 so the 
tensor mode is even deeper inside its respective horizon 
than the scalar mode is. The only concern is a penalty 
to computational efficiency as the modes become highly 
oscillatory when they deep within their horizon. 

With all the integration constants fixed, the full set of 
differential equations and p| can be inte¬ 

grated until both the scalar and tensor modes are well 
outside the sound and light horizons respectively. This 
requirement can be parametrised by a constant B <C 1. 
Following the same argument, if k = B aH , we have 
c s k < B aH as c s < 1. In summary we integrate the 
mode equations from a time such that c s k = AaH until 
k = B aH with A^> 1 and B 1. When calculating the 


bispectrum (for isosceles triangles) we have a third hori¬ 
zon to consider for the squeezed/folded mode. Similar 
arguments can be made and one should take care to en¬ 
sure all relevant modes exit their horizons and satisfy the 
relevant initial conditions, we have We found the bispec¬ 
trum to converge when A ~ 400 and B ~ 1/100. Higher 
values of A significantly increased the computation time 
due to the oscillatory nature of the mode functions while 
providing no real benefit. Smaller values of B did not 
affect the accuracy or computation time. 

With the scalar and tensor power spectra in hand, the 
observables n s and r can be calculated directly following 
their definitions, either as a function of scale k or at a 
specific “pivot” scale fc* for comparison with conventional 
models 


n s (k+) 


r{K) 


din (fc 3 P c ( fc )) 


din k 


k=k* 


Pd**) ’ 


( 20 ) 

( 21 ) 


where the factor of 8 in the definition of r arise from 
the definition of the tensor perturbations and from the 
fact that two independent polarisations contribute to the 
total power. 


B. Computation of the bispectrum 

The bispectrum of ( is the simplest, lowest-order mo¬ 
ment, where we expect to see deviations from a pure 
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FIG. 3: Dependence of /nl on S on shape and sound speed for the smallest scale k = & max = 10 -2 (Mpc) -1 . The 
optimum delta occurs when the relevant curve has converged. The left panel shows how the c) dependence varies for 
each shape evaluated at c~ 2 = 3. There is a mild shape dependence in the optimum S where squeezed triangles 
require smaller S values. As a consequence, noise from folded configurations occurs at larger values of S so the 
optimum S must lie between these two cases. The right panel shows how the S dependence varies with sound speed 
dependence evaluated in the equilateral limit. There is remarkably little dependence on c s even at very small sound 

speeds. 


Gaussian statistic. It corresponds to a tree-level three- 
point vertex for an interacting quantum field and will 
be the most dominant form of non-Gaussianity as higher 
order moments are expected to be suppressed by higher 
order terms in both the HSR parameters and level of cur¬ 
vature perturbations with ~ 10 -5 . In the isotropic 
limit it reduces to a function of three variables, the mag¬ 
nitudes of the wavevectors ki, k 2 , and k 3 making up the 
allowed, closed triangles in Fourier space 


(CkiCfcCks) = (27T) 3 5( 3 )(k 1 +k 2 + k 3 )B(fc 1 ,fc 2 ,fc3), ( 22 ) 

where the delta function imposes the closed triangle con¬ 
dition due to isotropy. We define the reduced, dimen¬ 
sionless, scale and shape dependent bispectrum as 


/nl(&i ,k 2 ,k 3 ) = 5B(k 1 ,k 2 ,k 3 )/ 
6(|Cfc 1 | 2 Kk a | 2 + |Cfc 1 | 2 Kk3| 2 + lCfc a | 2 |Cfc,| 2 ) , (23) 


This is different to the usual /nl, scale free, amplitude 
for the bispectrum quoted in the literature [28]. 


The weighting introduced in the definition of /nl (23) 
is known as the “local” weighting. Other definitions 
are used in the literature depending on the expected 
shape dependence of the signal. When observational con¬ 
straints are obtained from data, such as with Planck [5] 
the various choices of weighting are used to define lim¬ 
its on different types of /nl- These include equilateral 
and orthogonal weightings. The limits reported in [5] are 


ai = 0.8 ± 5.0 


fZ il = ~4 ±43. 


/ortho 

/NL 


-26 ± 21 . 


The most dominant contribution to the bispectrum 
comes from 0 expanded to third order in /. Follow¬ 
ing [31] [32] the third-order action for single field inflation 
with a constant sound speed c s is 


= M 2 / d 4 x 


/■ 


2 a 3 e 

3flcf 


72 " 1 ) C 3 + -5- ( A2- h 3 ( 1 “ 72 ) )CC 2 + -^(e+! - Cg)C(<9C) 2 


39 


ae, 2 a 6 e 


c 

3^3 


+ ^(.V-e)( d C 


(1 - |) CdiCdid 2 C + ^aPcdid- 2 Cdid- 2 c 


(24) 

(25) 


Section III.B of [T| discussed why the action is written in 


the form (24) in order to deal with apparent divergences 
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FIG. 4: Shape dependence of /nl for several trajectories evaluated at c~ 2 = 1 (left) and c~ 2 = 3 (right). All values 
of /nl are normalised to their value at (3 = 1 . Single field inflation models generically peak in the equilateral limit 
but because they must follow the consistency relation in the limit /3 0 their (3 dependence is much sharper. 


and we refer the reader to that work for further detail. A 
c s 7 ^ 1 , and indeed, an arbitrary time-dependant c s pro¬ 
vides no further complications in dealing with the third- 
order action. 

The ”In-In formalism” j3QH32] is used to calculate the 
bispectrum and ultimately /nl- Using ( [24] ) to define an 
interaction Hamiltonian and treating ChaO as a scalar 
field with canonical commutation relations, the bispec¬ 
trum can be reduced to a single integral over N. 



Here 5s[z\ denotes the imaginary part of the imaginary 
number z. Nq and N± represent times when the largest 
and smallest scales are sufficiently deep inside and far 
outside the sound-horizon respectively, using the same A 
and B parameters as described above. Z(N) implicitly 
depends on the shape and scale of the triangle but the 
function arguments have been omitted for brevity. 

We now specialise to the case where Aq = = k and 

ks = /3k where 0 < (3 < 2. This simple parametrisation 
covers many cases of interest. The squeezed, equilateral, 
and folded limits correspond to (3 = 0, 1 and 2 respec¬ 
tively. Z(N) then takes on the following form: 

Z(N) = (fiC 2 Cp + f2(% + fsCCCp + M%), 

fi =4 u, 

h = < 2+,?2) (si) (“+y , '- 3<) ) • < 2? > 

h = 12 « - ^ + (1 - /? 2 )e + (A - 1^ Aj , 

/ 4 = 6m - 2. (477 + 2(/3 2 - l)e + - 1^ ^ 2 e 2 ^ , 

where C = C k, C = d(/dN, (p = (p k and u = 1 - c~ 2 . 


At early times in the limit A —>> oo, \Z(N)\ —>> oo. How¬ 
ever we deform the integration contour by a small, imag¬ 
inary component iS so that the oscillations arising from 
(15) become exponentially suppressed. This is the usual 


choice of contour one makes when calculating interacting 
correlation functions. In this limit (15) becomes 


Cfc -4 lim 


^ ^ s ^—ic s k(l-\-i8)r 

s^o 2M p ia V ke 


(28) 


as r —)> — oo and the integral converges at very early 
times. 


1. Regulating the integral 


To calculate the bispectrum we integrate (26) numer¬ 
ically. Analytically, after performing the integral, one 
could take the limit 5 0 to obtain an answer that is 

well behaved. Unfortunately this is not possible numeri¬ 
cally and gives rise to large errors. We cannot integrate 
over an infinite range in time, i.e. from A = oo, a(N) = 0 
or N = — oo, so there will always be a sharp integration 
cutoff at very early times. Because of this sharp cutoff, 
the oscillations in the integrand result in large fluctua¬ 
tions in the final answer even though they should cancel 
out if the integration constant is formally extended to 


A solution o this problem is to add an exponential 
damping factor similarly to the one introduced in ( [28]) . 
This was the first approach taken by Chen et. al. in ]33| . 
However there are some issues with this method. Firstly 
the amplitude of the integrals tend to be suppressed re¬ 
sulting in an underestimation of the bispectrum. In ad¬ 
dition, the optimal value for the damping factor S needs 
to be fine tuned for each scale considered [33] . 

An alternative method exists which which does not suf¬ 
fer from these issues. It was first used in [34 and then 
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FIG. 5: c s dependence of /nl for several trajectories evaluated at (3 = 1 (left) and (3 = 0.1 (right). All values of /nl 
are normalised to their value at /3 = 1. For equilateral triangles the /3 dependence is much stronger. In the squeezed 
limit the c s dependence becomes much smaller but remains non-zero. 


expanded on in [1 . We refer the reader to [I] for the 
details. The method splits the integral into two parts at 
an arbitrary split point defined by c s k = XaH. X needs 
to be large enough for (15) to be a good approximation 
for all three modes. Some integration by parts is per¬ 
formed then X is chosen to minimise the error on the 
bispectrum. Unfortunately this method does not work 
for c s 7 ^ 1 because of the new £ /3 term. The method 
still prevents the contributions from the oscillations at 
early times from diverging but the £ /3 term still intro¬ 
duces a large oscillatory signature to the final integral. 
We therefore adopted the first method employing an im¬ 
proved exponential damping factor 


Z(N) ->■ Z(N)e~ 6 ^T , 


(29) 


in the numerical integration. 

Fig. [l] shows the dependence of /nl on the suppression 
factor S for n = 1 in both the squeezed and folded limits. 
For this figure, and all the other S dependence figures, 
a random trajectory was taken with s = 1.5, x = 1 and 
Lmax = 4, as defined in ( [l0| ). We see that if 5 is too 
small, the early time oscillations are not sufficiently sup¬ 
pressed producing a large amount of noise. This noise is 
exaggerated for large values of k. Secondly, if S is too 
large, the damping factor will interfere with the time de¬ 
pendence around the time of horizon crossing. This is 
the most dominant contribution to the integral so it will 
no longer be a good approximation to the bispectrum. 
For this choice of n = 1 it is hard to justify an optimal 
value of S where /nl has converged. 

Another issue is that the optimal S depends on the 
shape of the triangle. Indeed, between the folded and 
squeezed cases the optimal S drops by an order of mag¬ 
nitude. This dependence can be reduced by adjusting 
the value of n. c s k/aH is very large at early times and 
of order 1 during horizon crossing. Therefore increasing 
n will give stronger weighting to the damping factor at 


early times, while interfering less with the horizon cross¬ 
ing time. We found n = 3 to give the best results. The 
S dependence for n = 3 is shown in Fig. [2] 

Most of the residual noise arises from large k modes, 
particularly in the folded configuration. In contrast most 
results calculated in the equilateral configuration are rel¬ 
atively clean. Fig. [3] shows how the S dependence varies 
with shape factor /3 and c s in the equilateral configura¬ 
tion. Fig. [2] motivates a choice of S ~ 0.005 and we use 
this suppression factor along with n = 3 for the remain¬ 
der of our calculations. 


IV. RESULTS 

One way to test our numerical results for robustness 
and consistency is by comparison with the the squeezed 
limit consistency relation [SOI EH- For any single field 
inflation model the following limit must hold 

lim (CfciCfesCfes) (27r) 3 (5 (3) (ki +k 2 +k 3 ) (n s - 1) P kl P k3 

k3<tiki,k2 

(30) 


or in our notation 

lim /nl K - 1) . (31) 

/3 — ^0 12 

It is important to emphasise here that this holds for 
all single field models independent of the value of c s or 
the prior we choose for the initial conditions of the back¬ 
ground trajectories. However increasing the value of c s 
or the HSR parameters typically increases the amplitude 
of /nl therefore we don’t necessarily expect all models to 
tend to the squeezed limit at the same rate. For example 
/3 = 0.1 might be “squeezed enough” for low values of c s 
but not for higher values. We first analyse the shape and 
sound speed dependence of the trajectories, elaborating 
















on the consistency relation in section [TV C| Unless stated 
otherwise, the trajectories are taken from a prior with 
x = 1, s = 1.5 and L max = 4. 


A. Shape dependence 

Fig. @ compares the shape dependence of trajectories 
evaluated at k m [ n = 10 -5 (Mpc ) -1 normalised to their 
equilateral values. As expected, for trajectories with 
shape dependence |/nl| peaks in the equilateral config¬ 
uration. As c s reduces, the amplitude of |/nl| typically 
increases but the trajectories must still obey the squeezed 
limit consistency relation where |/nl| ^ 10 -2 . This ex¬ 
aggerates the shape dependence of all the trajectories, 
even those which appear fiat when c s = 1 . 

It is worth noting that in the squeezed limit, the shape 
dependence is curved in comparison to the roughly lin¬ 
ear dependence in the folded limit. This is in agreement 
with m where the authors show that corrections linear 
in /? drop out. Any terms linear in k 3 must contract sym¬ 
metrically with the remaining two modes. As they have 
equal magnitudes in opposite directions they will can¬ 
cel out leaving only quadratic corrections in k%. In the 
folded limit this cancellation does not occur producing 
the linear dependence shown in Fig. [4| 


B. c s dependence 


Fig. [5] compares the dependence of /nl on c s for equi¬ 
lateral and squeezed triangles. These values are nor¬ 
malised to their values at c s = 1. To a good approxima¬ 
tion the dependence is linear in c ~ 2 and much stronger for 
equilateral triangles. This shows that for fixed /3 = 0.1 
one can still obtain large /nl by choosing an arbitrarily 
small c s . At c s = 1 , /nl is typically small and negative 
so as c s —* 0 /nl becomes large and positive. The close 
linear dependence on c ~ 2 is not surprising and it clearly 
arises from the functions fi in (27). 


C. Monte Carlo Plots 

The scale dependence is linear to a good approximation 
and can easily be analysed. To this end we define n/ NL 
as 


^/nl (^*) P') 


d/NL(&,/3) 


dlnfc 


k —/c* 


(32) 


As discussed in [45H471 it is possible to define a scale 
dependence as long as the shape of the triangle is kept 
fixed. Our definition is different to the usual definition 
of n / NL which is the derivative of | In /nl | and this is 
simply to avoid difficulties arising when /nl ~ 0. Recall, 
reducing c s often in induces a sign change as can be seen 
in Fig. [5] 


Fig-© shows numerous Monte Carlo plots for various 
sound speeds and shapes. Each plot consists of 2 18 tra¬ 
jectories with their colour representing n/ NL . The top 
two figures show that all trajectories tend towards the 
squeezed limit consistency relation even for small sounds 
speeds c s < 1. The consistency relation 5 (n s — 1)/12 
is shown by the red-dashed line. To reach the consis¬ 
tency relation in the c~ 2 = 3 case, a much smaller (3 
was required (and consequently the value of 5 had to be 
lowered, recall Fig. [ 2 ]). 

In the equilateral case one can see clearly how a small 
sound speed deforms the inflationary attractor. For ex¬ 
ample in the c s = 1 case, the consistency relation acts 
as a firm upper limit for /nl- The deviation from the 
consistency relation is simply proportional to e > 0 and 
f(k) > 0 defined in [30] . A small c s < 1 clearly violates 
this relation deforming the distribution significantly, re¬ 
sulting in a large positive /nl- In the folded limit, the 
distribution is reduced back again to be parallel with the 
consistency relation, although this time with a positive, c s 
dependent offset. 

To illustrate the flexibility of the method Fig. [7] shows 
a distribution with c~ 2 = 100 with colour of the tra¬ 
jectories now representing the third slow roll parameter 
£ = 2 A evaluated shortly after horizon crossing and the 
tensor-to-scalar ratio r. The dashed lines represent the 
current Planck constraints on n s = 0.968 db 0.006 [4]. 
Planck also constrains f^ 11 = — 4 ± 43 [5] although it is 
important to remember that there is not an exact one-to- 
one correspondence between our /nl calculated here and 
the one constrained by Planck [5] due to assumptions on 
scale-invariance. 

For example in power law inflation £ = e 2 and is often 
assumed to be vanishingly small. However at these sound 
speeds, one can see that a small variation in £ can lead 
to an appreciable change in /nl even though it is likely 
to be neglected. 

From the right panel in Fig. [7] one can also see that 
for small c s tighter constraints on r require larger |/nl|- 
From one perspective this is not surprising as, to lead¬ 
ing order, r ~ 16c s e m so smaller sounds speeds nat¬ 
urally induce smaller r. However one has to remember 
that the right panel in Fig. [7] shows trajectories for fixed 
c s = 0.1. The changes in /nl and r can only be induced 
by the slow-roll parameters (and H). More concretely 
smaller values of e are thus expected to produce more 
non-Gaussianity. This is in contrast to the c s = 1 case 
where larger values in e produce more non-Gaussianity. 
Indeed it is often quoted that /nl 00 e. From the plots 
this is fairly easy to explain. Increasing e always con¬ 
tributes negatively to /nl- It just so happens that at 
c s = 1 , /nl is small and negative so they add construc¬ 
tively. On the other hand reducing c s always contributes 
positively to /nl eventually inducing a sign change. As 
soon as /nl changes sign, increasing e reduces the amount 
of non-Gaussianity. 
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V. DISCUSSION 

We have outlined a full, numerical calculation of the 
bispectrum with a particular emphasis on single field 
models of inflation with non-canonical speed of sound. 
The calculation is challenging due to the oscillatory na¬ 
ture of the integrands involved which is exacerbated for 
the case with c s ^ 1 and we have shown how regularising 
the integrals can lead to stable results with the correct 
choice of numerical damping terms. The methods ex¬ 
plored in this work can be used to investigate the scale 
and shape dependence of the bispectrum signal produced 
by an epoch of inflation. 

For convenience we have adopted a more general de¬ 
scription of bispectrum signal than that normally quoted 
in the literature by re-defining a scale and shape depen¬ 
dent /nl? which always tends to 5 (n s — 1)/12 as the shape 
parameter (3 0. For lower values of c s , |/nl| is typi¬ 

cally much greater and thus requires much smaller values 
of (3 to recover the squeezed limit consistency relation. 

If future observational surveys of the CMB or large 
scale structure become accurate enough to constrain any 


scale dependence of the non-Gaussian signal then our 
work could be applied to the calculation of accurate 
model of the bispectrum to be used in likelihood eval¬ 
uations of the data. This is not currently possible as 
the strongest limit on non-Gaussianity come from an 
ad-hoc analysis of Planck CMB maps assuming a scale- 
independent and fixed shape templates for the bispec¬ 
trum leading to constraints on a single amplitude pa¬ 
rameter. Whilst these results may be consistent with the 
simplest model of inflation, if a non-zero amplitude for 
/nl were ever to be measured, more accurate parametri- 
sations of the non-Gaussianity will be useful to try to 
gain a better understanding of the nature of the infla- 
ton and its connection with extensions to the standard 
model of particle physics. This will particularly become 
a priority if primordial tensor modes are not discovered 
at levels r ~ 0.01 — 0.1. 
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FIG. 6: Monte Carlo plots for c s = 1 (left) and c~ 2 = 3 (right). From top to bottom the shape configurations are 
evaluated in the squeezed, equilateral and folded limits respectively. The red-dashed line represents the consistency 
relation 5 (n s — 1)/12. The colour of each trajectory illustrates the scale dependence of the bispectrum, nj NL . For 
squeezed c~ 2 = 3 (top-right) it was necessary to reduce p = 0.02 to recover the squeezed limit as opposed to j3 = 0.1 
for the c s = 1 case. This increased computation time by roughly an order of magnitude. 
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FIG. 7: Monte Carlo plot for a very small sound speed c~ 2 = 100 evaluated in the equilateral limit. The red-dashed 
lines represent the recent Planck constraints on n s [4]. The left panel shows how small variations in £ can change 
/nl- The right panel shows the corresponding tensor-to-scalar ratio r for each trajectory. 



















